import pathlib
import sys

# make coastpy library importable by appending root to path
cwd = pathlib.Path().resolve()
proj_dir = cwd.parent.parent
sys.path.append(str(proj_dir / "src"))

# definitly not a best practice, but a workaround to avoid many warnings caused by shapely 2.0 release
import io
import socket
import warnings

import cartopy.crs as ccrs
import geopandas as gpd
import geoviews as gv
import geoviews.tile_sources as gts
import geoviews.tile_sources as gvts
import holoviews as hv
import hvplot.pandas
import hvplot.xarray
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import panel as pn
import requests
import rioxarray
import shapely
import xarray as xr
from tqdm import tqdm

from coastpy.geometries import geo_bbox

warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)

/var/folders/p8/qxnqn2ns5kbczrp91kx4_1nm0000gn/T/ipykernel_17541/1733808755.py:7: UserWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:

import os
os.environ['USE_PYGEOS'] = '0'
import geopandas

In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
  import geopandas as gpd

DATA_DIR = proj_dir / "data"
coastsat_metadata_pacific_fp = DATA_DIR / "03_CoastSat_metadata_layer_Pacific.geojson"
metadata = gpd.read_file(coastsat_metadata_pacific_fp)

for c in ["Tier1Images", "Tier2Images", "S2Images", "pk"]:
    metadata[c] = metadata[c].astype("i8")

for c in ["TideRange", "WaveHeight", "RelativeTideRange"]:
    metadata[c] = metadata[c].astype("f8")

metadata.head()

	Tier1Images	Tier2Images	S2Images	TideRange	WaveHeight	RelativeTideRange	pk	geometry
0	410	113	207	1.65	1.22	1.35	0	POINT (180.00000 -16.15293)
1	256	33	214	1.79	0.60	2.97	1	POINT (179.16992 -16.41845)
2	255	56	415	1.86	0.81	2.29	2	POINT (178.51402 -16.78987)
3	301	78	214	1.47	0.60	2.43	3	POINT (179.18860 -16.72214)
4	495	176	214	1.51	1.01	1.49	4	POINT (179.89177 -16.65793)

Part 1: compute relative tidal range with respect to waves

todo: check with Kilian which waves heights he used mean_tidal_range / mean_wave_height

metadata["RelativeTideRange"] = metadata["TideRange"] / metadata["WaveHeight"]
# metadata["RelativeTideRange"] = ...  # your code to compute tidal range over the wave height

pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Tide versus wave dominated coasts",
        style={"color": "black"},
        width=500,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

w_threshold = pn.widgets.FloatSlider(
    name="Threshold", start=0.1, end=10.0, step=0.1, value=1
)


@pn.depends(w_threshold.param.value)
def plot(threshold):
    wave_dominated = metadata[metadata["RelativeTideRange"] < threshold].copy()
    tide_dominated = metadata[metadata["RelativeTideRange"] > threshold].copy()

    return (
        gts.EsriImagery()
        * wave_dominated.hvplot(geo=True, color="blue", label="Wave dominated")
        * tide_dominated.hvplot(geo=True, color="green", label="Tide dominated")
    ).opts(legend_position="bottom_right")


app = pn.Column(
    title_bar,
    pn.Row(w_threshold),
    pn.Row(plot),
)
app

from ipyleaflet import Map, basemaps

m = Map(basemap=basemaps.Esri.WorldImagery, scroll_wheel_zoom=True)
m.center = 37.768137, -122.511066
m.zoom = 16
m.layout.height = "800px"
m

bbox = geo_bbox(m.west, m.south, m.east, m.north)

bbox.explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

transects = gpd.read_file(DATA_DIR / "03_CoastSat_transect_layer.geojson")
transects_roi = gpd.sjoin(transects, bbox)

transects_roi.explore()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[10], line 1
----> 1 transects_roi.explore()

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/geodataframe.py:1981, in GeoDataFrame.explore(self, *args, **kwargs)
   1978 @doc(_explore)
   1979 def explore(self, *args, **kwargs):
   1980     """Interactive map based on folium/leaflet.js"""
-> 1981     return _explore(self, *args, **kwargs)

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/explore.py:366, in _explore(df, column, cmap, color, m, tiles, attr, tooltip, popup, highlight, categorical, legend, scheme, k, vmin, vmax, width, height, categories, classification_kwds, control_scale, marker_type, marker_kwds, style_kwds, highlight_kwds, missing_kwds, tooltip_kwds, popup_kwds, legend_kwds, map_kwds, **kwargs)
    363         map_kwds["max_zoom"] = tiles.get("max_zoom", 18)
    364         tiles = tiles.build_url(scale_factor="{r}")
--> 366 m = folium.Map(
    367     location=location,
    368     control_scale=control_scale,
    369     tiles=tiles,
    370     attr=attr,
    371     width=width,
    372     height=height,
    373     **map_kwds,
    374 )
    376 # fit bounds to get a proper zoom level
    377 if fit:

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/folium/folium.py:274, in Map.__init__(self, location, width, height, left, top, position, tiles, attr, min_zoom, max_zoom, zoom_start, min_lat, max_lat, min_lon, max_lon, max_bounds, crs, control_scale, prefer_canvas, no_touch, disable_3d, png_enabled, zoom_control, **kwargs)
    272     zoom_start = 1
    273 else:
--> 274     self.location = validate_location(location)
    276 Figure().add_child(self)
    278 # Map Size Parameters.

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/folium/utilities.py:78, in validate_location(location)
     71         raise ValueError(
     72             "Location should consist of two numerical values, "
     73             "but {!r} of type {} is not convertible to float.".format(
     74                 coord, type(coord)
     75             )
     76         )
     77     if math.isnan(float(coord)):
---> 78         raise ValueError("Location values cannot contain NaNs.")
     79 return [float(x) for x in coords]

ValueError: Location values cannot contain NaNs.

TIMEOUT_SEC = 5  # default timeount in seconds

# Get the transect id's by saving the transects in our region of interest as a list.
transect_ids = transects_roi.TransectId.to_list()

# Here we download the data per transect, but if it takes too long we default to Ocean beach, California.
try:
    data = []
    for transect_id in tqdm(transect_ids):
        url = url = f"http://coastsat.wrl.unsw.edu.au/time-series/{transect_id}/"
        # Pandas.read_csv() can handle url's but it doesn't accept a timeout argument,
        # so we download the data via requests.
        r = requests.get(url, timeout=TIMEOUT_SEC).content
        s = pd.read_csv(
            io.StringIO(r.decode("utf-8")),
            header=None,
            names=["date", "shoreline_position"],
        )
        s["date"] = pd.to_datetime(s["date"])
        s = s.set_index("date").replace({"None": np.nan}).astype("f8")
        s.name = transect_id
        data.append(s)
    column_names = [s.name for s in data]
    shorelines = pd.concat(data, axis=1)
    shorelines.columns = column_names

except:
    # If it takes too long to get the data from the server, default to Ocean Beach, California.
    # First set the map to that area, so that we have the correct bbox.
    m.center = 37.768137, -122.511066
    m.zoom = 16
    m.layout.height = "800px"
    bbox = geo_bbox(m.west, m.south, m.east, m.north)
    transects_roi = gpd.sjoin(transects, bbox)
    shorelines_fp = DATA_DIR / "03_shorelines_ocean_beach.parquet"
    shorelines = pd.read_parquet(shorelines_fp)

  0%|                                         | 0/8 [00:00<?, ?it/s]

shorelines.head()

	usa_CA_0215-0003	usa_CA_0215-0004	usa_CA_0215-0005	usa_CA_0215-0006	usa_CA_0215-0007	usa_CA_0215-0008	usa_CA_0215-0009	usa_CA_0215-0010	usa_CA_0215-0011	usa_CA_0215-0012	usa_CA_0215-0013	usa_CA_0215-0014	usa_CA_0215-0015	usa_CA_0215-0016	usa_CA_0215-0017	usa_CA_0215-0018
date
1984-05-02 18:13:15	71.4	87.6	103.9	109.0	115.4	NaN	120.1	124.8	NaN	NaN	NaN	NaN	148.7	146.1	131.7	NaN
1984-06-03 18:14:10	74.4	70.1	98.5	116.8	117.9	141.0	127.4	134.9	NaN	166.6	NaN	NaN	NaN	139.3	133.6	NaN
1984-06-10 18:20:29	62.9	71.4	90.4	118.9	111.7	137.4	140.6	126.9	137.5	158.1	150.3	150.1	155.8	150.8	142.7	NaN
1984-06-26 18:20:42	60.7	74.3	91.6	112.9	127.7	133.8	152.4	156.0	157.7	171.2	164.9	174.4	165.9	148.2	144.7	NaN
1984-09-07 18:16:09	77.3	91.1	110.3	128.9	141.6	152.9	166.0	173.3	173.5	173.8	174.8	184.2	195.9	202.3	190.7	185.3

check this for inspiration on how to clor the plot> https://foundations.projectpythia.org/core/xarray/enso-xarray.html

pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Historical shoreline series",
        style={"color": "black"},
        width=500,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

options = transects_roi.TransectId.to_list()
transect_slider = pn.widgets.Select(
    name="Transect", options=options, value=np.random.choice(options)
)


@pn.depends(transect_slider.param.value)
def plot_transects(transect_id):
    transect = transects_roi.loc[transects_roi["TransectId"] == transect_id].copy()

    return (
        gts.EsriImagery()
        * transects_roi.hvplot(geo=True, color="blue")
        * transect.hvplot(geo=True, color="red")
    )


@pn.depends(transect_slider.param.value)
def plot_series(transect_id):
    shoreline_selected = shorelines.hvplot(x="date", y=[transect_id], color="red")
    shorelines_ = shorelines[
        shorelines.columns.difference([transect_slider.value])
    ].hvplot(x="date", alpha=0.4)

    return shorelines_ * shoreline_selected


app = pn.Column(
    title_bar,
    pn.Row(transect_slider),
    pn.Row(pn.Column(plot_transects), plot_series),
)
app

mei_fp = DATA_DIR / "03_meiv2.data"
mei = (
    pd.read_csv(mei_fp, delim_whitespace=True, index_col=0, header=None)
    .reset_index()
    .rename({0: "year"}, axis="columns")
    .melt(id_vars="year", var_name="month")
)
mei["date"] = pd.to_datetime(mei["year"].astype(str) + "-" + mei["month"].astype(str))
mei = (
    mei.drop(["year", "month"], axis="columns")
    .set_index("date")
    .sort_index()
    .replace({-999: np.nan})
)
mei.hvplot(x="date")

mei_interp = mei.reindex(shorelines.index, method="nearest")

shorelines_ = shorelines.copy()
shorelines_.index = mei_interp.value.to_list()
shorelines_.hvplot(kind="scatter")

shorelines_ = shorelines.diff().copy()
shorelines_.index = mei_interp.value.to_list()
shorelines_.hvplot(kind="scatter", alpha=0.2, color="blue")